﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary
{
    /// <summary>
    /// This class provides methods to interpolate a classical cubic spline out of specified points.
    /// </summary>
    [Serializable]
    public class CubicSplineInterpolation
    {
        /// <summary>
        /// The points for the interpolation.
        /// </summary>
        private readonly Point2DCollection interpolationPoints;

        /// <summary>
        /// Initializes a new instance of the <see cref="CubicSplineInterpolation"/> class.
        /// </summary>
        /// <param name="interpolationPoints">The points for the interpolation.</param>
        public CubicSplineInterpolation(Point2DCollection interpolationPoints)
        {
            if (interpolationPoints == (Point2DCollection) null)
            {
                throw new ArgumentNullException("interpolationPoints");
            }

            this.interpolationPoints = interpolationPoints;
        }

        /// <summary>
        /// Interpolates the cubic splines with the specified points.
        /// </summary>
        /// <returns>The interpolated cubic spline.</returns>
        public CubicSpline Interpolate()
        {
            if (this.interpolationPoints.Count < 4)
            {
                throw new ArgumentException("Minimum four points were needed for an interpolation.");
            }

            CubicSpline result = new CubicSpline();
            double[] h = new double[this.interpolationPoints.Count - 1];
            double[] y = new double[this.interpolationPoints.Count - 2];

            this.interpolationPoints.Sort();

            for (int i = 1; i < this.interpolationPoints.Count; i++)
            {
                h[i - 1] = this.interpolationPoints[i].X - this.interpolationPoints[i - 1].X;
            }

            for (int i = 2; i < this.interpolationPoints.Count; i++)
            {
                double a = (this.interpolationPoints[i].Y - this.interpolationPoints[i - 1].Y) / h[i - 1];
                double b = (this.interpolationPoints[i - 1].Y - this.interpolationPoints[i - 2].Y) / h[i - 2];

                y[i - 2] = 6 * (a - b);
            }

            Matrix solveMatrix = new Matrix(y.Length, y.Length + 1);

            solveMatrix.SetValueAtPosition(new MatrixPosition(0, 0), 2 * (h[0] + h[1]));
            solveMatrix.SetValueAtPosition(new MatrixPosition(0, 1), h[1]);
            solveMatrix.SetValueAtPosition(new MatrixPosition(0, y.Length), y[0]);

            for (int i = 1; i < y.Length - 1; i++)
            {
                solveMatrix.SetValueAtPosition(new MatrixPosition(i, i - 1), h[i]);
                solveMatrix.SetValueAtPosition(new MatrixPosition(i, i), 2 * (h[i] + h[i + 1]));
                solveMatrix.SetValueAtPosition(new MatrixPosition(i, i + 1), h[i + 1]);
                solveMatrix.SetValueAtPosition(new MatrixPosition(i, y.Length), y[i]);
            }

            solveMatrix.SetValueAtPosition(new MatrixPosition(y.Length - 1, y.Length - 1), 2 * (h[0] + h[1]));
            solveMatrix.SetValueAtPosition(new MatrixPosition(y.Length - 1, y.Length - 2), h[1]);
            solveMatrix.SetValueAtPosition(new MatrixPosition(y.Length - 1, y.Length), y[y.Length - 1]);

            double[] tempuri = solveMatrix.SolveLinearSystem();
            double[] solveResults = new double[tempuri.Length + 2];

            solveResults[0] = 0;
            solveResults[solveResults.Length - 1] = 0;

            for (int i = 1; i < solveResults.Length - 1; i++)
            {
                solveResults[i] = tempuri[i - 1];
            }

            double[] resultsD = this.CalculateComponentD();
            double[] resultsB = CalculateComponentB(solveResults);
            double[] resultsA = CalculateComponentA(solveResults, h);
            double[] resultsC = this.CalculateComponentC(solveResults, h);
            SimplePolynomial[] polynomialsA = new SimplePolynomial[resultsD.Length];
            SimplePolynomial[] polynomialsB = new SimplePolynomial[resultsD.Length];
            SimplePolynomial[] polynomialsC = new SimplePolynomial[resultsD.Length];
            SimplePolynomial[] polynomialsD = new SimplePolynomial[resultsD.Length];
            LinearFactorCollection[] collectionA = new LinearFactorCollection[resultsD.Length];
            LinearFactorCollection[] collectionB = new LinearFactorCollection[resultsD.Length];

            for (int i = 0; i < resultsD.Length; i++)
            {
                collectionA[i] = new LinearFactorCollection();
                collectionB[i] = new LinearFactorCollection();
            }

            for (int i = 0; i < resultsD.Length; i++)
            {
                for (int j = 0; j < 3; j++)
                {
                    collectionA[i].Add(new LinearFactor(-this.interpolationPoints[i].X));
                }

                for (int l = 0; l < 2; l++)
                {
                    collectionB[i].Add(new LinearFactor(-this.interpolationPoints[i].X));
                }

                polynomialsC[i] = new SimplePolynomial(-this.interpolationPoints[i].X, 1);
            }

            for (int i = 0; i < polynomialsA.Length; i++)
            {
                polynomialsA[i] = collectionA[i].MultiplyFactors();
                polynomialsB[i] = collectionB[i].MultiplyFactors();
            }

            for (int i = 0; i < polynomialsA.Length; i++)
            {
                polynomialsA[i].MultiplyWithScalar(resultsA[i]);
                polynomialsB[i].MultiplyWithScalar(resultsB[i]);
                polynomialsC[i].MultiplyWithScalar(resultsC[i]);
            }

            SimplePolynomial[] endResults = new SimplePolynomial[resultsD.Length];

            for (int i = 0; i < resultsD.Length; i++)
            {
                endResults[i] = polynomialsA[i] + polynomialsB[i] + polynomialsC[i]
                                + new SimplePolynomial(resultsD[i]);
            }

            for (int i = 0; i < resultsD.Length; i++)
            {
                result.Add(new CubicSplineComponent(endResults[i].Coefficients[0],
                                                    endResults[i].Coefficients[1],
                                                    endResults[i].Coefficients[2],
                                                    endResults[i].Coefficients[3],
                                                    this.interpolationPoints[i], this.interpolationPoints[i + 1]));
            }

            return result;
        }

        /// <summary>
        /// Calculates the component D of the polynomial.
        /// </summary>
        /// <returns>The calculated components.</returns>
        private double[] CalculateComponentD()
        {
            double[] result = new double[this.interpolationPoints.Count - 1];

            for (int i = 0; i < result.Length; i++)
            {
                result[i] = this.interpolationPoints[i].Y;
            }

            return result;
        }

        /// <summary>
        /// Calculates the component B of the polynomial.
        /// </summary>
        /// <param name="s">The s components.</param>
        /// <returns>The calculated components.</returns>
        private static double[] CalculateComponentB(double[] s)
        {
            double[] result = new double[s.Length - 1];

            for (int i = 0; i < result.Length; i++)
            {
                result[i] = s[i] / 2;
            }

            return result;
        }

        /// <summary>
        /// Calculates the component A of the polynomial.
        /// </summary>
        /// <param name="s">The s components.</param>
        /// <param name="h">The h components.</param>
        /// <returns>The calculated components.</returns>
        private static double[] CalculateComponentA(double[] s, double[] h)
        {
            double[] result = new double[s.Length - 1];

            for (int i = 0; i < result.Length; i++)
            {
                result[i] = (s[i + 1] - s[i]) / (6 * (h[i]));
            }

            return result;
        }

        /// <summary>
        /// Calculates the component C of the polynomial.
        /// </summary>
        /// <param name="s">The s components.</param>
        /// <param name="h">The h components.</param>
        /// <returns>The calculated components.</returns>
        private double[] CalculateComponentC(double[] s, double[] h)
        {
            double[] result = new double[s.Length - 1];

            for (int i = 0; i < result.Length; i++)
            {
                result[i] = ((this.interpolationPoints[i + 1].Y - this.interpolationPoints[i].Y) / h[i]) -
                            (((2 * h[i] * s[i]) + (h[i] * s[i + 1])) / 6);
            }

            return result;
        }
    }
}